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ABSTRACT 

To  reduce  the  influence  of  atmospheric  turbulence  on  images  of  space-based  objects  we  are  developing  a  maximum 
a  posteriori  deconvolution  approach.  In  contrast  to  techniques  found  in  the  literature,  we  are  focusing  on  the 
statistics  of  the  point-spread  function  (PSF)  instead  of  the  object.  We  incorporated  statistical  information  about  the 
PSF  into  multi-frame  blind  deconvolution.  Theoretical  constraints  on  the  average  PSF  shape  come  from  the  work  of 
D.  L.  Fried  while  for  the  univariate  speckle  statistics  we  rely  on  the  gamma  distribution  adopted  from  radar/laser 
speckle  studies  of  J.  W.  Goodman.  Our  aim  is  to  develop  deconvolution  strategy  which  is  reference-less,  i.e.,  no 
calibration  PSF  is  required,  extendable  to  longer  exposures,  and  applicable  to  imaging  with  adaptive  optics.  The 
theory  and  resulting  deconvolution  framework  were  validated  using  simulations  and  real  data  from  the  3.5m 
telescope  at  the  Starfire  Optical  Range  (SOR)  in  New  Mexico. 

1.  INTRODUCTION 


Solutions  to  recovery  of  high-resolution  images  when  observing  through  atmospheric  turbulence  usually  fall  into  the 
software  (“post -processing”)  or  the  hardware  (adaptive  optics,  interferometry)  category  or  the  combination  of  both. 
Even  with  very  expensive  adaptive  optics  (AO)  systems  it  is  necessary  to  use  deconvolution  (image  reconstruction) 
to  remove  image  blurring  completely  [1,2], 

Successful  restoration  of  images  degraded  by  atmospheric  turbulence  was  first  achieved  using  the  so-called  speckle 
imaging  techniques  [3],  Speckle  imaging  works  because  average  short -exposure  power  spectrum  has  non-negligible 
spectral  content  up  to  the  telescope’s  diffraction  limit  D/I  [4],  On  the  other  hand,  for  long  exposures  average  optical 
transfer  function  (OTF)  quickly  drops  to  values  below  the  noise  limit  above  the  cut-off  rfl.  No  information  can  be 
recovered  from  the  part  of  the  spectrum  where  signal-to-noise  ratio  (SNR)  falls  below  unity.  It  is  this  SNR  cutoff 
that  actually  determines  available  resolution,  with  or  without  deconvolution,  and  the  SNR  cutoff  is  nearly  always 
short  of  the  diffraction  cutoff  [5],  Therefore,  deconvolution  of  long -exposure  images  taken  without  AO  cannot  result 
in  reliable  amplification  of  the  high-frequency  content. 

Multi-frame  blind  deconvolution  (MFBD)  [6-8]  is  an  image  reconstruction  method  relying  on  the  availability  of 
several  images  of  an  object.  In  addition,  many  of  the  MFBD  algorithms  rely  on  short  exposures  for  the  reason  stated 
above.  The  multiplicity  of  image  frames  acts  as  an  implied  constraint  because  the  object  is  common  to  every  image. 
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while  noise  and  turbulence  fluctuations  vary  randomly  between  frames  [5],  Even  with  this  advantage,  MFBD  can 
get  easily  trapped  in  local  minima  [9], 

Regularization  is  the  most  popular  approach  to  balancing  resolution  enhancement  and  amplification  of  noise.  Early 
termination  of  the  iterative  process  or  an  additional  spectral  filter  can  be  thought  of  as  simple  regularization  methods 
[5],  In  our  work  we  rely  on  the  Bayesian  formulation  of  the  problem  of  finding  a  true  object  o  and  a  (stochastic)  PSF 
h  which  together  with  noise  generated  the  recorded  data  i: 


p(o,hli) 


p{\  lo,h)-p(o)-p(h) 

P( i) 


(i) 


The  theorem  states  that  the  conditional  probability  of  the  object  being  equal  to  o  and  the  PSF  being  equal  to  h  given 
that  we  recorded  data  i  is  equal  to  the  product  of:  conditional  probability  of  the  data  taking  on  the  value(s)  i  given  o 
and  h,  probability  of  the  object  being  equal  to  o,  and  the  probability  of  the  PSF  being  equal  to  h,  divided  by  the 
probability  of  obtaining  the  data  i  which  is  always  taken  to  be  unity.  In  the  maximum  a  posteriori  (MAP)  framework 

one  finds  estimates  for  the  object,  6  ,  and  the  PSF(s),  h  ,  which  jointly  maximize  /?(<>,  h  i): 

[6, h]  =  argmax p{ o, h|  i)  =  argmax  p(i\ o, h)  x  p( o)  x  p(h)  (2) 

o,  h  o, h 


It  is  often  useful  to  rewrite  the  above  equation  in  terms  of  the  negative  log -likelihoods.  Then,  6  can  be  defined  as 
the  object  that  minimizes  a  compound  criterion  Joh( i): 


JoM  (i)  =  Ji  (o,  h)  +  Ja  (o)  +  J  h  (h) 


(3) 


where  the  negative  log-likelihoods  are  defined  according  to  the  rule  Jx  (z)  =  -ln/?(x,y|  z).  Since  a  probability 

density  function  p(o)  for  the  object  is  in  general  not  known,  solutions  in  the  form  of  functions  with  desirable 
mathematical  properties  (e.g.  noise  suppression,  edge  enhancement)  are  often  used.  As  these  ad  hoc  formulas  are  not 
real  probability  density  functions  (PDFs),  regularization  parameters  must  be  used  to  balance  their  influence  on  the 
cost  function  in  Equation  (2).  Often  these  parameters  have  to  be  chosen  manually  [1,8].  Apart  from  the  problem  of 
choosing  the  right  value  for  them,  the  mere  form  of  the  prior  (e.g.  an  object’s  power  spectral  density)  could  be 
applicable  only  to  a  limited  class  of  real  objects,  although  the  object’s  power  spectral  density  could  be  estimated 
from  the  images  themselves  [10], 

For  these  reasons  we  focus  on  PSF  statistics  p(h).  This  part  of  Equation  (2)  was  almost  always  removed  from  the 
MAP  approach  [7,10].  Formulas  for  wavefront  statistics  and,  by  extension,  for  short-exposure  PSF  statistics  have 
been  used  in  the  field  of  wavefront  reconstruction  and  deconvolution  from  wavefront  sensing  [11]  but  not  in  the 
field  of  image  restoration.  The  few  papers  which  deal  with  the  PSF  prior  do  so  for  the  case  of  very  long  exposures 
[1,12],  Here  we  present  ideas  on  how  to  incorporate  the  body  of  knowledge  about  turbulence -induced  statistics  of 
intensity  into  the  MAP  framework. 

2.  CONSTRAINTS  ON  THE  TURBULENT  PSF 


There  exists  a  great  deal  of  theoretical  information  about  the  turbulent  PSF  that  has  not  yet  been  used  in  multi -frame 
blind  deconvolution.  For  example,  average  ensemble  PSF  is  completely  specified  by  known  optical  parameters,  such 
as  observing  wavelength  and  the  aperture  of  the  sensor,  and  one  unknown  parameter  describing  the  integrated  effect 
of  the  atmospheric  turbulence  between  the  source  and  the  observer.  The  often-used  parameterization  of  this 
integrated  effect  is  through  the  Fried’s  coherence  length  r0  [13],  If  r0  at  the  time  of  the  observations  could  be 
measured  then  it  could  be  used  to  generate  an  average  PSF  for  subsequent  deconvolution.  As  mentioned  in  the 
introduction,  deconvolution  of  the  average  long-exposure  image,  even  with  the  perfectly  known  PSF,  cannot  result 


in  a  diffraction-limited  image  when  SNR  is  below  unity  for  high  spatial  frequencies.  Unfortunately,  this  is  very  often 
the  case  for  faint  space  objects. 

On  the  other  hand,  short  exposures  allow  the  recovery  of  the  diffraction-limited  image.  Additionally,  a  sequence  of 
short  exposures  of  an  arbitrary  object  can  be  processed  to  yield  r0.  In  order  to  remove  the  object  being  viewed  from 
the  image  formation  equation  an  object-cancelling  transformation  is  used  [14].  Recently,  we  proposed  an  original 
transformation,  which  we  call  “Fourier  contrast’’  method  [15].  Here  we  give  a  brief  description  how  it  works. 

Image  formation  equation,  expressed  in  the  Fourier  domain,  is: 

I(w)  =  0(l7)H(i7)  (4) 


where  u  is  a  spatial  frequency  vector  in  the  Fourier  plane  and  I(;7) ,  O (u)  and  H(;7)  stand  for  Fourier  transforms 
of  the  (instantaneous)  image,  object,  and  speckle  PSF,  respectively.  For  each  frequency  we  calculate  the  mean 
values  and  standard  deviations  of  the  power  spectra: 


c;(«) 


v  ar(|l(f7)|2 ) 1/2 

H*> 


0(t7)|2  var(|H(t7)|2)1/2  _  var(|H(t7)|2)l/2 
|0(w)|2  (|H(t7)|2)  (|H(7)|2) 


(5) 


where  (}j  denotes  average  and  var(.)  denotes  variance.  We  use  the  term  “contrast”,  and  denote  it  with  letter  C  as  is 
customary  in  research  pertaining  to  speckle  [4],  As  can  be  appreciated  from  Equation  (5)  in  the  absence  of  noise  the 
object  disappears  in  C;  (;7)  .  With  noise  present  in  the  data  the  object  does  not  cancel  out  after  the  transformation 
but  its  influence  is  very  small  on  the  low-frequency  part  of  C7  (u)  which  is  the  part  where  rQ  estimation  must  be 

performed.  Models  for  ^|H(t7)|  ''j  and  var(|H(t7)|  )have  been  developed  using  the  theory  of  partially-developed 

speckle  [15].  Since  the  models  depend  only  on  D,  X  and  r0  it  is  possible  to  estimate  Fried’s  parameter  by  fitting 
theoretical  templates  to  measured  C;  (u)  . 

Once  r0  is  given  the  long-,  and  short-exposure  (tip-tilt  corrected)  PSFs  can  be  generated  using  the  theory  of  Fried 
[13],  In  this  theory,  the  optical  transfer  function  OTF,  i.e.,  the  Fourier  transform  of  the  PSF,  for  the  two  cases  of 
long  exposures  and  registered  short  exposures  is  given  by: 


H0  {it )  x  H  u:  ( u ) 

(6) 

Hj  (w  )  —  H0  {tt  )x  Hs,  {tt ) 

(7) 

where  H  u  (it )  represents  the  average  long -exposure  OTF  of  the  atmosphere,  H  SE  (u )  is  the  average  short-exposure 
OTF  of  the  atmosphere,  H ;  (it )  and  H s,  (it)  are  the  overall  long-  and  short-exposure  average  OTFs  (including  the 
effect  of  the  telescope  diffraction).  For  a  diffraction -limited  circular  aperture  of  diameter  D  we  have: 


H0(m)  =  — 
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f  i _ 7  i . 
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where  u  =  u  ,  A  is  the  average  wavelength,  and  z  is  the  distance  from  the  exit  pupil  to  the  image  plane. 


Fried  developed  expressions  for  H/7  (</ )  and  H  SE  (u ) : 


exd  -3.44 
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where  a  is  a  parameter  that  varies  between  1/2  when  there  are  both  intensity  and  phase  variations  across  the 
collecting  aperture  and  1  when  only  phase  distortions  are  present.  By  Fourier-transforming  Equations  (6)  or  (7)  one 
obtains  the  average  PSF  (h(x,  yjj  . 

To  solve  Equation  (2)  we  need  to  know  the  probability  density  function  of  turbulent  PSF  p(h).  In  this  paper  we 
present  the  equations  for  the  statistics  of  intensity,  what  we  call  the  “focal-plane”  constraint,  but  we  will  also 
mention  at  this  stage  that  it  is  possible  to  constrain  the  statistics  of  the  PSF  in  pupil  plane  and  in  the  Fourier  plane 
[16]. 

With  infinitesimally  small  detectors  (pixels),  exposure  times  and  filter  bandwidths  one  can  model  p(h)  as  the 
exponential  distribution  [4]: 


/;(")= firm 


(h(x,y)) 


h  (k,x,y) 
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where  K  is  the  number  of  speckle  frames  and  X  and  Y  are  the  numbers  of  pixels  in  x  and  y  directions.  Speckle 
patterns  in  general  will  be  integrated:  spatially,  temporally,  and  spectrally.  To  model  these  effects,  use  is  made  of  the 
gamma  PDF  [4]: 
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where  M  is  the  general  integration  parameter,  a  product  of  the  spatial,  spectral  and  temporal  integration  parameters 

[17]. 

In  closing  this  section  we  also  give  an  equation  for  the  data-fidelity  term  in  Equation  (2).  Under  the  assumption  of 
spatially-stationary  Gaussian  noise  with  standard  deviation  o  we  have: 
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(i|o,h)=(« 
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(13) 


We  assume  p(o)  is  a  uniform  distribution  and  solve  the  MAP  problem  using  Equations  (12)  and  (13). 


3.  DISTRIBUTION  TESTING 


Because  real  data  contains  both  Poisson  and  readout  components  on  top  of  speckle  variability,  the  assumption  of  the 
gamma  model  for  the  integrated  speckle  was  tested  against  a  noise-free  simulation.  The  simulation  employs  only 
one  phase  screen  and  assumes  isoplanatic  conditions  (only  global  image  motion).  Scintillation  effects  are  ignored. 

Independent  phase  screens  are  generated  using  the  classic  FFT -based  method  [18],  whereby  an  array  of  random 
numbers  is  filtered  according  to  the  von-Karman  spectrum  and  then  inverse  Fourier -transformed.  The  outer  scale 
was  set  to  25  m.  We  deliberately  did  not  use  subharmonics  correction  and  we  removed  global  tip/tilts  from  the 
wavefronts  in  order  to  have  the  same  -  in  the  statistical  sense  -  speckles  falling  on  a  given  pixel.  Size  of  the 
telescope  aperture  was  set  to  50  cm  and  256x256  pixel  arrays  were  used  for  phase  screens  giving  pupil  sampling  of 
2  mm  per  pixel.  Fried’s  parameter  was  set  to  2  cm  at  500  nm  which  is  also  the  reference  wavelength  of  the 
simulations.  The  phase  screens  were  converted  to  phasors,  multiplied  by  a  circular  aperture  and  embedded  in  an 
array  of  zeros  of  size  512x512  pixels.  The  field  was  then  Fourier-transformed.  Square  of  the  modulus  of  the  result 
gives  the  PSF.  Pixel  scale  of  the  images  is  49  Lirad  (Nyquist  sampling  at  the  shortest  wavelength).  The  simulation  is 
polychromatic:  images  corresponding  to  ten  wavelengths  between  500  and  700  nm  are  generated  assuming  linear 
wavefront  scaling.  The  resulting  ten  PSFs  are  co-added  and  the  result  constitutes  one  polychromatic  noise-free  PSF. 
Additionally,  five  such  PSFs  were  added  to  simulate  temporal  averaging.  This  binning  size  was  chosen  by  trial-and- 
error  with  the  goal  to  obtain  similar  values  of  integration  parameter  M  across  the  field-of-view  as  in  the  SOR  data.  In 
the  end,  a  sequence  of  200  integrated  images  was  generated  out  of  1000  instantaneous  polychromatic  frames. 

The  Kolmogorov-Smirnov  [19,20]  and  Anderson-Darling  [21]  statistical  tests  were  used  to  obtain  quantitative 
confidence  levels  for  the  null  hypothesis  that  a  speckle  sample  from  a  given  pixel  comes  from  the  gamma 
distribution.  In  both  tests  the  test  statistic  D  is  a  measure  of  distance  between  the  empirical  distribution  function 
(EDF)  -  which  is  the  proportion  of  the  observed  values  that  are  less  than  or  equal  to  a  particular  value  -  and  the 
hypothesized  cumulative  density  function  (CDF).  In  the  Kolmogorov-Smirnov  test  D  is  the  maximum  of  the 
absolute  difference  between  EDF  and  CDF,  while  in  the  Anderson-Darling  test  D  corresponds  to  the  weighted 
squared  difference  between  EDF  and  CDF.  If  the  null  hypothesis  -  that  CDF  is  the  underlying  distribution  -  is 
correct,  EDF  should  be  close  to  CDF  and  D  should  be  close  to  zero.  The  result  of  the  test,  the  p- value,  gives  the 
probability  that  a  value  of  D  at  least  as  large  as  the  one  observed  would  have  occurred  if  the  null  hypothesis  were 
true.  The  greater  the  p- value  the  more  confidence  one  can  have  in  the  null  hypothesis.  Usually,  values  higher  than 
0.05  are  taken  as  strong  evidence  that  a  sample  came  from  the  hypothesized  distribution  (p- value  is  bounded 
between  0  and  1). 

When  all  the  parameters  of  a  hypothesized  CDF  are  specified  a  priori,  there  exist  approximate  formulas  for  the 
computation  of  the  p- value.  Unfortunately  this  was  not  the  case  here:  both  (h(.r,  yjj  and  M  in  Equation  (12)  have  to 

be  estimated  from  the  sample.  Because  of  that  the  bootstrap  simulation  had  to  be  used  in  order  to  obtain  the  p- values 
[22],  The  parameters  of  the  gamma  distribution  were  estimated  using  the  method  of  moments.  In  this  approach, 
sample  moments  are  equated  to  the  unobservable  population  moments.  Then  the  equations  relating  the  distribution 
parameters  to  the  population  moments  are  solved.  Specifically,  computation  of  M  follows  as: 


(h(*»  y)) 

<h(*,  y)) 


Computation  of  (h(.r,  yj)  is  trivial.  With  both  these  parameters  gamma  CDF  can  be  computed  and  the 
Kolmogorov-Smirnov  and  Anderson-Darling  tests  can  be  applied  to  the  sample. 

The  bootstrap  simulation  is  a  Monte-Carlo  approach  to  estimating  the  confidence  intervals  for  the  null  hypothesis. 
The  random  variables  with  the  hypothesized  distribution  are  generated,  and  the  parameters  of  the  distribution  are 
estimated  from  the  simulated  samples  using  the  same  method  as  for  the  observed  sample.  Again,  the  test  statistic  D 
is  calculated.  The  result  of  the  bootstrap  method  is  the  number  of  times  the  test  statistic  D  calculated  from  the 
generated  sample  is  greater  than  or  equal  to  D  calculated  from  the  observed  sample.  This  number  divided  by  the 
total  number  of  simulations  gives  the  p- value.  Gamma-distributed  random  variables  were  generated  using  the 
method  of  Marsaglia  and  Tsang  [23].  Ten  thousand  samples  of  the  same  size  as  a  given  data  set  were  generated.  The 


p- values  were  averaged  for  all  pixels  having  the  same  distance  from  the  center  of  the  PSF.  This  is  because  the 
spectral  component  of  parameter  M  increases  with  distance  from  the  center  (the  spatial  and  temporal  components  are 
constant).  One  can  explain  it  in  the  following  way:  M  can  be  thought  of  as  a  number  of  speckles  contributing  to  one 
pixel.  Imagining  a  discrete  set  of  increasing  wavelengths,  producing  radially  expanding  PSFs,  we  see  that  stepping 
over  wavelengths  results  in  new  speckles  traversing  a  pixel  under  observation  if  it  is  away  from  the  PSF  center.  The 
further  a  pixel  is  from  the  center  the  more  speckles  will  traverse  it  (compare  the  sharp  speckles  close  in  to  the 
smooth  streaks  further  out  from  the  center  in  the  left  panel  of  Figure  2). 

The  results  of  statistical  testing  are  shown  in  Figure  1  for  the  distances  of  0-30  pixels  from  the  center.  In  almost  all 
cases  we  obtained  values  higher  than  0.3.  This  indicates  that  the  gamma  PDF  provides  a  good  model  for  the 
integrated  speckle  (spatially,  temporally  and  spectrally). 


0.5 


0.4 


0.3 


*  6 
o  a  ^  a 


.  o 

A  A  O  o  O 
A  a  A  <? 


0aA<>oaa  o 
aa8  *  <>A  °° 


A  ° 


0.2  - 


0.1  - 


O  Kolmogorov— Smirnov  test 
A  Anderson -Darling  test 


0.0 


5  10 

Distance  from  the  image  center  (\/D) 


15 


Fig.  1.  Confidence  levels,  quantified  by  the  K-S  or  A-D  p-value,  that  the  integrated  speckle  intensity  is  governed  by  the  gamma  PDF. 


4.  RESULTS 


The  model  PDFs  from  Equations  (12)  and  (13)  were  converted  to  negative  log-likelihoods  and,  together  with 
analytic  gradients,  inserted  into  a  Variable  Metric  with  Limited  Memory  (VMLM)  optimizer  [24],  I-band  (800-900 
nm)  observations  of  the  bright  single  star  HR2219  (Figure  2,  top  left  panel)  have  been  obtained  with  the  3.5  m 
telescope  at  the  Starfire  Optical  Range  with  adaptive  optics  switched  off  (but  with  the  tip/tilt  system  switched  on). 
Thousand  frames  were  recorded,  with  exposure  time  set  to  10  ms.  We  obtained  D/r0  =  34  from  the  Fourier  contrast 
method. 

To  execute  our  MFBD  on  real  data  we  needed  to  know  the  integration  parameter  M  to  be  plugged  into  Equation 
(12).  The  estimation  of  M  for  spatial-only  integration  is  straightforward  when  one  knows  the  auto-correlation 
function  of  speckle  and  pixel  shape  [4,17].  Also,  we  conjecture  that  estimation  of  spectral  M  could  be  done  easily 
under  the  assumption  of  rectangular  filter  and  the  similarity  of  spatial  and  spectral  auto -correlation  functions  (the 
latter  will  be  the  scaled  version  of  the  former  with  the  scaling  parameter  depending  on  the  distance  from  the  image 
center;  see  the  text  preceding  Figure  1).  Both  these  components  of  the  total  M  would  have  to  be  computed  only  once 
for  a  given  telescope-camera-filter  combination.  The  problematic  part  is  the  estimation  of  temporal  M.  Use  will  have 
to  be  made  of  the  analytical  temporal  auto-correlation  function  which  has  a  dependency  on  the  vertical  profile  of 
wind  velocity  [25]. 

Before  we  attempt  a  fully  reference-less,  analytical  computation  of  spatial-spectral-temporal  M,  we  estimate  it  from 
the  combination  of  real  data  statistics  and  corresponding  noise-less  simulation.  Figure  2  shows  M  (in  the  right 
panels)  of  either  real  speckle  images  (top  left  panel)  or  the  simulated  ones  (bottom  left  panel).  Simulated  data 


corresponding  to  SOR  parameters  (wavelength,  filter  bandwidth,  pixel  size,  telescope  diameter,  etc.)  was  generated 
using  the  approach  described  in  Section  3,  but  with  subharmonics  correction  of  the  phase-screens.  Computation  of 
M  at  each  location  in  the  PSF  was  done  by  the  means  of  Equation  (14).  Note  how  the  M  image  for  real  speckles  is 
only  valid  where  the  SNR  is  high,  i.e.,  where  the  speckle  fluctuations  are  higher  than  noise.  We  used  the  M  image  of 
simulated  speckles  to  scale  M  of  SOR  data  by  a  constant  value  to  correct  for  the  difference  in  the  integration  time 
between  the  simulated  (instantaneous  exposure)  and  the  real  (10  ms)  images.  In  such  way  we  have  been  able  to  use 
the  gamma  prior  for  the  whole  field-of-view  which  was  not  possible  in  the  case  of  the  real  M  image  only  (Figure  2, 
top  right  panel)  because  of  noise. 


Fig  2.  Top  left:  a  real  speckle  image  taking  with  the  3.5  m  SOR  telescope,  exposure  time  equal  10  ms,  spectral  bandwidth  from  800  to  900  nm, 
pixel  size  80.5  nrad.  Top  right:  the  associated  M  image.  Bottom  left:  a  simulated  speckle  image  for  the  same  telescope,  instantaneous  exposure 
time,  bandwidth  from  800  to  900  nm,  pixel  size  equal  80.5  nrad.  Bottom  right:  the  corresponding  M  image. 


Speckle  images  from  the  single  star  HR2219  were  used  as  PSFs.  They  were  convolved  with  a  schematic 
representation  of  the  Hubble  Space  Telescope  (HST,  Figure  3)  which  acted  as  the  true  object.  The  resulting  images 
were  corrupted  with  Gaussian  noise  of  standard  deviations  equal  to  1  or  5  counts. 


Fig  3.  True  object  used  for  the  simulations. 

Multi  frame  blind  deconvolution  with  gamma  PSF  prior  was  used  to  deconvolve  the  simulated  data.  Results  are 
shown  in  Figure  4.  Reconstructions  show  significantly  enhanced  contrast  with  respect  to  either  a  single  frame  or  to 
the  shift-and-add  image.  Figure  5  shows  evolution  with  the  number  of  iterations  of  the  peak  signal-to-noise  ratio 
(PSNR)  defined  as: 


L2 


PSNR  =  101og10 


MSE 


(15) 


where  L  is  the  dynamic  range  of  the  image.  PSNR  is  calculated  using  the  mean  squared  error  (MSE)  which 
compares  a  particular  reconstruction  A  with  the  reference  true  image  B: 


MSE = X  Z  ( A(* ’  y)  -  >?))2 

Xl  X=i  y=l 


(16) 


From  Figure  5  it  is  clear  that  use  of  the  gamma  distribution  prior  on  the  PSF  prevents  noise  amplification  in  the 
reconstruction. 

Finally,  in  Figure  6  we  compare  the  estimated  PSFs  with  the  real  ones  used  for  the  simulations.  Our  estimates  follow 
the  general  morphology  and  shape  of  the  real  PSFs  in  terms  of  size,  elongation  and  maximum-peak  positions  but 
lack  the  same  high-frequency  content. 


Fig  4.  Top  row:  case  of  readout  noise  with  a  =  1 .  Bottom  row:  case  of  readout  noise  with  a  =  5.  Left  column:  single  frames.  Middle  column: 
shift-and-add  images.  Right  column:  final  reconstructions.  Linear  scale  from  0  to  250  is  the  same  in  all  panels. 


Fig  5.  PSNR  values  vs.  number  of  iterations  for  two  noise  levels. 


Fig  6.  Top  row:  real  speckle  images  used  as  PSFs.  Bottom  row:  Estimated  PSFs  from  MFBD. 


5.  OUTLOOK 

Currently,  we  are  estimating  the  integration  parameter  M  using  the  PSFs  and  simulations.  To  satisfy  the  “no 
reference”  goal  of  the  project  we  will  have  to  estimate  M  analytically  with  the  assumed  spatial,  spectral  and 
temporal  auto-correlations.  The  validity  of  these  assumptions  will  have  to  be  tested  against  real  data.  Temporal  and 
spectral  integration  effects  are  particularly  interesting.  Modeling  of  them  has  been  attempted  in  the  context  of 
speckle  interferometry  but  the  resulting  formulas  are  rather  complex  and  have  dependencies  on  several  parameters 
which  cannot  be  easily  obtained  from  the  data,  like  phase  variance  [26].  We  hope  to  provide  an  alternative,  simple 
description  of  the  effect  of  integration  on  various  speckle  statistics.  Subsequently,  we  will  test  PSF  priors  expressed 
in  the  pupil  plane  and  in  the  Fourier  plane  [16]. 


ACKNOWLEDGEMENTS 

Effort  sponsored  by  U.S.  Air  Force  Office  of  Scientific  Research,  Air  Force  Material  Command,  under  grant 
FA8655-12-1-21 15,  and  the  Spanish  Ministry  of  Science  under  grant  AyA2008-01225.  The  U.S.  Government  is 
authorized  to  reproduce  and  distribute  reprints  for  Governmental  purpose  notwithstanding  any  copyright  notation 
thereon. 


REFERENCES 

[1]  L.  M.  Mugnier,  T.  Fusco,  and  J.-M.  Conan,  "MISTRAL:  a  myopic  edge -preserving  image  restoration  method, 
with  application  to  astronomical  adaptive -optics-corrected  long -exposure  images,"  J.  Opt.  Soc.  Am.  A  21, 
1841-1854(2004). 

[2]  R.  Baena  Galle,  J.  Nunez,  and  S.  Gladysz,  "Extended  object  reconstruction  in  adaptive -optics  imaging:  the 
multiresolution  approach,"  Astronomy  &  Astrophysics,  555,  A69-A84,  (2013).. 

[3]  J.  C.  Dainty,  "Stellar  speckle  interferometry",  in  Laser  Speckle  and  Related  Phenomena,  edited  by  J.  C.  Dainty, 
Springer  Verlag,  Heidelberg,  2nd  Edition  (1984). 

[4]  J.  W.  Goodman,  "Speckle  phenomena  in  optics:  theory  and  applications,"  (Roberts  &  Company,  2006). 


[5]  D.  W.  Tyler,  "Deconvolution  of  adaptive  optics  image  data,"  Proceedings  of  the  Center  for  Adaptive  Optics, 

(2001). 

[6]  T.  J.  Schulz,  "Multiframe  blind  deconvolution  of  astronomical  images,"  J.  Opt.  Soc.  Am.  A  10,  1064-1073 
(1993). 

[7]  D.  G.  Sheppard,  B.  R.  Hunt,  and  M.  W.  Marcellin,  "Iterative  multiframe  superresolution  algorithms  for 
atmospheric-turbulence-degraded  imagery,"  J.  Opt.  Soc.  Am.  A  15,  978-992  (1998). 

[8]  C.  L.  Matson,  K.  Borelli,  S.  Jefferies,  C.  C.  Beckner,  Jr.,  E.  K.  Hege,  and  M.  Lloyd-Hart,  "Fast  and  optimal 
multiframe  blind  deconvolution  algorithm  for  high-resolution  ground-based  imaging  of  space  objects,"  Appl. 
Opt.  48,  A75-A92  (2009). 

[9]  J.  Nagy  and  V.  Mejia-Bustamante,  "MFBD  and  the  local  minimum  trap,"  Proceedings  of  the  AMOS 
Conference,  2009,  Ed.:  S.  Ryan,  The  Maui  Economic  Development  Board,  p.ElO 

[10]  L.  Blanco  and  L.  M.  Mugnier,  "Marginal  blind  deconvolution  of  adaptive  optics  retinal  images,"  Opt.  Express 
19,23227-23239  (2011). 

[11]  L.  M.  Mugnier,  C.  Robert,  J.-M.  Conan,  V.  Michau,  and  S.  Salem,  "Myopic  deconvolution  from  wave-front 
sensing,"  J.  Opt.  Soc.  Am.  A  18,  862-872  (2001). 

[12]  A.  MacDonald,  S.  C.  Cain,  and  E.  E.  Armstrong,  "Maximum  a  posteriori  image  and  seeing  condition 
estimation  from  partially  coherent  two-dimensional  light  detection  and  ranging  images,"  Opt.  Eng.  45,  086201- 
1-13  (2006). 

[13]  D.  L.  Fried,  "Optical  resolution  through  a  randomly  inhomogeneous  medium  for  very  long  and  very  short 
exposures,"  J.  Opt.  Soc.  Am.  56,  1372-1379  (1966). 

[14]  O.  von  der  Liihe,  "Estimating  Fried’s  parameter  from  a  time  series  of  an  arbitrary  resolved  object  imaged 
through  atmospheric  turbulence,"  J.  Opt.  Soc.  Am.  A  1,  510-519  (1984). 

[15]  S.  Gladysz,  R.  Baena  Galle,  R.  Johnson,  and  L.  Kann,  "Image  reconstruction  of  extended  objects: 
demonstration  with  the  Starfire  Optical  Range  3.5m  telescope,"  Proceedings  of  SPIE,  Volume  8535,  85350M- 
1-13(2012). 

[16]  R.  Baena  Galle,  S.  Gladysz,  L.  Mugnier,  V.  Gudimetla,  R.  Johnson,  and  L.  Kann,  "Physically-constrained 
Multi-frame  Blind  Deconvolution,"  in  Imaging  and  Applied  Optics,  J.  Christou  and  D.  Miller,  eds.,  OSA 
Technical  Digest  (online)  (Optical  Society  of  America,  2013),  paper  JW2A.3 

[17]  S.  E.  Skipetrov,  J.  Peuser,  R.  Cerbino,  P.  Zakharov,  B.  Weber,  and  F.  Scheffold,  "Noise  in  laser  speckle 
correlation  and  imaging  techniques,"  Opt.  Express  18,  14519-14534  (2010). 

[18]  B.  L.  McGlamery,  "Computer  simulation  studies  of  compensation  of  turbulence  degraded  images,"  in  Image 
Processing,  J.  C.  Urbach,  ed„  Proc.  SPIE  74,  225-233  (1976). 

[19]  A.  N.  Kolmogorov,  "Sulla  determinazione  empirica  di  una  legge  di  distribuzione,"  Giornale  dell'  Istituto 
Italiano  degli  Attuari  4,  83-91  (1933). 

[20]  N.  V.  Smirnov,  "Table  for  estimating  the  goodness  of  fit  of  empirical  distributions,"  The  Annals  of 
Mathematical  Statistics,  19,  279-281  (1948). 

[21]  T.  W.  Anderson  and  D.  A.  Darling,  "Asymptotic  theory  of  certain  "goodness-of-fit"  criteria  based  on  stochastic 
processes,"  Annals  of  Mathematical  Statistics  23,  193-212  (1952). 

[22]  S.  Gladysz,  J.  C.  Christou,  L.  W.  Bradford,  and  L.  C.  Roberts,  Jr.,  "Temporal  variability  and  statistics  of  the 
Strehl  ratio  in  adaptive-optics  images,"  PASP  120,1132-1143  (2008). 

[23]  G.  Marsaglia  and  W.  Tsang,  "A  simple  method  for  generating  gamma  variables,"  ACM  Transactions  on 
Mathematical  Software  26,  363-372  (2000). 

[24]  E.  Thiebaut,  "Optimization  issues  in  blind  deconvolution  algorithms,"  in  Astronomical  Data  Analysis  II,  SPIE 
Vol.  4847,  174-183  (2002). 

[25]  F  Roddier,  J  M  Gilli  and  G  Lund,  "On  the  origin  of  speckle  boiling  and  its  effects  in  stellar  speckle 
interferometry,"  J.  Optics  13,  263-271  (1982). 

[26]  J.  Ohtsubo,  "Effects  of  finite  spectral  bandwidth  and  focusing  error  on  the  transfer  function  in  stellar  speckle 
interferometry,"  J.  Opt.  Soc.  Am.  A  2,  667-673  (1985). 


